home *** CD-ROM | disk | FTP | other *** search
/ USGS: Oil & Gas Potential…National Wildlife Refuge / USGS - Oil & Gas Potential of the Arctic National Wildlife Refuge - Disc 2.iso / mac / MEcode / MESamp.for < prev    next >
Text File  |  1999-02-11  |  9KB  |  733 lines
  1. c   program name: MESamp.for.
  2.  
  3. c    Tests user supplied pairwise correlation structure
  4.  
  5. c     for internal consistency, performs Cholesky decomposition &
  6.  
  7. c     generates correlated sample numbers based upon rank correlation.
  8.  
  9. c     The resultant sample numbers will be used to aggregate ANWR play
  10.  
  11. c     results.
  12.  
  13. c
  14.  
  15. c    Written 4/2/98, Schuenemeyer
  16.  
  17. c
  18.  
  19. c    Input file:
  20.  
  21. c      Unit 10 - DepA.dat, Contains play number, number of plays, play 
  22.  
  23. c           names and correlation matrix
  24.  
  25. c    Output file:
  26.  
  27. c      Unit 11 - Rand.dat, contains sample index numbers
  28.  
  29. c
  30.  
  31. c    Subroutines required: jacobi.for dcholdc.for buble.for
  32.  
  33. c      & bublei.for (included)
  34.  
  35. c
  36.  
  37.       parameter(epsilon=0.001,lx=50)
  38.  
  39.       character pln*15,title*70
  40.  
  41.     real*8 wk2(50),dtri(50,50)
  42.  
  43.     real mult
  44.  
  45.       dimension rho(50,50),tri(50,50),nply(50)
  46.  
  47.      1 ,asvd(50,50),ev(50),v(50,50),wk1(lx),id(lx)
  48.  
  49.      2 ,iv(50),pln(50),u(50),au(50),x(10000,11)
  50.  
  51.     3 ,is(10000),st(10000),it(10000)
  52.  
  53.       write(*,2)
  54.  
  55.     2 format(' Program MESamp.for generates correlated sample indices')
  56.  
  57. c   open correlation matrix file
  58.  
  59.       open(10,file='DepA.dat',status='old')
  60.  
  61.     open(11,file='Rand.dat')
  62.  
  63.       do i=1,np
  64.  
  65.        id(i)=i
  66.  
  67.       end do
  68.  
  69.       nerr=0
  70.  
  71. c
  72.  
  73. c   read correlation file
  74.  
  75.     read(10,'(a70)') title
  76.  
  77. c      np is number of plays
  78.  
  79.     read(10,*) np        
  80.  
  81.     do i = 1,np
  82.  
  83. c      iv is play num, nply is num of replications, pln is play name
  84.  
  85.      read(10,'(i2,1x,i5,1x,a15)')iv(i),nply(i),pln(i)
  86.  
  87.     end do
  88.  
  89. c   read "correlation matrix"
  90.  
  91.     do i=1,np
  92.  
  93.        read(10,*) (rho(i,j),j=1,i)
  94.  
  95.     end do
  96.  
  97.     do i=1,np
  98.  
  99.        do j=1,i 
  100.  
  101.         rho(j,i)=rho(i,j)
  102.  
  103.         asvd(i,j)=rho(i,j)
  104.  
  105.         asvd(j,i)=rho(i,j)
  106.  
  107.         tri(i,j)=0.0
  108.  
  109.         tri(j,i)=0.0
  110.  
  111.      end do
  112.  
  113.       end do
  114.  
  115. c   compute eigenvalues; see if `corr matrix' pos semi-def
  116.  
  117.       do i=1,np
  118.  
  119.        rho(i,i)=1.0
  120.  
  121.        asvd(i,i)=1.0
  122.  
  123.        tri(i,i)=0.0
  124.  
  125.       end do
  126.  
  127.       call jacobi(asvd,np,lx,ev,v,nrot)
  128.  
  129.       do i=1,np
  130.  
  131.        wk1(i)=ev(i)
  132.  
  133.       end do
  134.  
  135.       call buble(wk1,id,np)
  136.  
  137.       write(*,15)
  138.  
  139.    15 format(/' Eigenvalues'/)
  140.  
  141.       do i=1,np
  142.  
  143.        write(*,17)i,wk1(i)
  144.  
  145.    17  format(i5,f7.3,5f7.3)
  146.  
  147.       end do
  148.  
  149. c   allows user to view eigenvalues
  150.  
  151.       pause 0
  152.  
  153.       if (wk1(1).ge.0.0) then
  154.  
  155.        write(*,19) wk1(1)
  156.  
  157.    19 format(/' Minimum eigenvalue=',f10.3/
  158.  
  159.      1 3x,'Matrix is correlation matrix-bias factor not needed')
  160.  
  161.      do i=1,np
  162.  
  163.       dtri(i,i)=1.0
  164.  
  165.       do j=1,i
  166.  
  167.        dtri(i,j)=rho(i,j)
  168.  
  169.        dtri(j,i)=dtri(i,j)
  170.  
  171.       end do
  172.  
  173.      end do
  174.  
  175.      goto 40
  176.  
  177.       end if
  178.  
  179.       bias=abs(wk1(1))+epsilon
  180.  
  181.       write(*,22)wk1(1),bias
  182.  
  183.    22 format(/' Min. eigenvalue=',f10.3,'. Bias=',f10.3)
  184.  
  185.       do i=1,np
  186.  
  187.        ev(i)=ev(i)+bias
  188.  
  189.       end do
  190.  
  191. c  Generate new adjusted correlation matrix
  192.  
  193.       do i=1,np
  194.  
  195.       do j=1,np
  196.  
  197.         do k=1,np
  198.  
  199.          tri(i,j)=tri(i,j)+v(i,k)*ev(k)*v(j,k)
  200.  
  201.         end do
  202.  
  203.       end do
  204.  
  205.       end do
  206.  
  207. c  Normalize new correlation matrix.
  208.  
  209.       xnc=tri(1,1)
  210.  
  211.       do i=2,np
  212.  
  213.        do j=1,i-1
  214.  
  215.         tri(i,j)=tri(i,j)/xnc
  216.  
  217.         tri(j,i)=tri(i,j)
  218.  
  219.      end do
  220.  
  221.       end do
  222.  
  223.     do i=1,np
  224.  
  225.     tri(i,i)=1.0
  226.  
  227.     dtri(i,i)=1.0
  228.  
  229.      do j=1,i
  230.  
  231.       dtri(i,j)=tri(i,j)
  232.  
  233.       dtri(j,i)=dtri(i,j)
  234.  
  235.      end do
  236.  
  237.     end do
  238.  
  239. c
  240.  
  241. c    Compute Cholesky decomposition to get triangular matrix
  242.  
  243.     call dcholdc(dtri,np,lx,wk2)
  244.  
  245.       write(*,18)
  246.  
  247.    18 format(/' New Correlation Matrix-upper right - '
  248.  
  249.      1 'Cholosky-lower left '/)
  250.  
  251.     do i=1,np
  252.  
  253.     write(*,33)iv(i),(dtri(i,j),j=1,np)
  254.  
  255.    33 format(i3,10f6.3)
  256.  
  257.       end do
  258.  
  259. c   Perform check on Cholesky decomp
  260.  
  261.    40 tdmax=0.
  262.  
  263.     do i=1,np
  264.  
  265.      do j=i,np
  266.  
  267.        t=0
  268.  
  269.        do k = 1,i
  270.  
  271.         t=t+dtri(j,k)*dtri(i,k)
  272.  
  273.        end do
  274.  
  275.        if(i.eq.j) then
  276.  
  277.        td=abs(t-1.)
  278.  
  279.        else
  280.  
  281.        td=abs(t-dtri(i,j))
  282.  
  283.        end if
  284.  
  285.        if(td.gt.tdmax) dtmax=td
  286.  
  287.      end do
  288.  
  289.     end do
  290.  
  291.     if (abs(td).gt.0.001) then
  292.  
  293.       write(*,*) ' Warning,Cholesky decomp suspect'
  294.  
  295.       write(*,*) '   Deviations too large', td
  296.  
  297.     end if
  298.  
  299. c       
  300.  
  301. c   Generate random correlated sample numbers
  302.  
  303. c
  304.  
  305.      iseed=871
  306.  
  307.     call seed(iseed)
  308.  
  309. c   num is number of samples
  310.  
  311.     num=10000
  312.  
  313.     do m=1,num
  314.  
  315.      do j=1,np
  316.  
  317. c   generates uniforms in range (-1,+1)
  318.  
  319.       u(j)=rand()*2. - 1.
  320.  
  321.        end do
  322.  
  323. c   multiply lower triangular * u to get dependent structure
  324.  
  325.        do i=1,np
  326.  
  327.       au(i)=0.0
  328.  
  329.       do j=1,i
  330.  
  331.        au(i)=au(i)+dtri(i,j)*u(j)
  332.  
  333.       end do
  334.  
  335.       x(m,i)=au(i)
  336.  
  337.        end do
  338.  
  339.       end do
  340.  
  341. c   rank the sample
  342.  
  343.     do i=1,np
  344.  
  345.      do m=1,num
  346.  
  347.       is(m)=m
  348.  
  349.       it(m)=m
  350.  
  351.       st(m)=x(m,i)
  352.  
  353.      end do
  354.  
  355.      call buble(st,is,num)
  356.  
  357.      call bublei(is,it,num)
  358.  
  359.      mult=1.0
  360.  
  361.      if(nply(i).gt.10000) mult=real(nply(i))/10000.
  362.  
  363.      do m=1,num
  364.  
  365.       x(m,i)=int(real(it(m))*mult)
  366.  
  367.      end do
  368.  
  369.     end do
  370.  
  371. c this is for Nig(many), play 11
  372.  
  373.       do m=1,num
  374.  
  375.      x(m,11)=int(x(m,10)*1.5432)
  376.  
  377.     end do
  378.  
  379. c write sample index numbers
  380.  
  381.     do m=1,num
  382.  
  383.      write(11,58)m,(x(m,j),j=1,np+1)
  384.  
  385.    58  format(i6,11f8.0)
  386.  
  387.       end do
  388.  
  389.       stop
  390.  
  391.       end
  392.  
  393.     SUBROUTINE jacobi(a,n,np,d,v,nrot)
  394.  
  395. C     Numerical Recipes Software, 1986-92
  396.  
  397. c      d is eigenvalaues of a in first n elements.
  398.  
  399. c      v contains the normalized eignevectors.
  400.  
  401. c      note d=v trans * a * v
  402.  
  403.       INTEGER n,np,nrot,NMAX
  404.  
  405.       REAL a(np,np),d(np),v(np,np)
  406.  
  407.       PARAMETER (NMAX=500)
  408.  
  409.       INTEGER i,ip,iq,j
  410.  
  411.       REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
  412.  
  413.       do 12 ip=1,n
  414.  
  415.         do 11 iq=1,n
  416.  
  417.           v(ip,iq)=0.
  418.  
  419. 11      continue
  420.  
  421.         v(ip,ip)=1.
  422.  
  423. 12    continue
  424.  
  425.       do 13 ip=1,n
  426.  
  427.         b(ip)=a(ip,ip)
  428.  
  429.         d(ip)=b(ip)
  430.  
  431.         z(ip)=0.
  432.  
  433. 13    continue
  434.  
  435.       nrot=0
  436.  
  437.       do 24 i=1,50
  438.  
  439.         sm=0.
  440.  
  441.         do 15 ip=1,n-1
  442.  
  443.           do 14 iq=ip+1,n
  444.  
  445.             sm=sm+abs(a(ip,iq))
  446.  
  447. 14        continue
  448.  
  449. 15      continue
  450.  
  451.         if(sm.eq.0.)return
  452.  
  453.         if(i.lt.4)then
  454.  
  455.           tresh=0.2*sm/n**2
  456.  
  457.         else
  458.  
  459.           tresh=0.
  460.  
  461.         endif
  462.  
  463.         do 22 ip=1,n-1
  464.  
  465.           do 21 iq=ip+1,n
  466.  
  467.             g=100.*abs(a(ip,iq))
  468.  
  469.             if((i.gt.4).and.(abs(d(ip))+
  470.  
  471.      1       g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then
  472.  
  473.               a(ip,iq)=0.
  474.  
  475.             else if(abs(a(ip,iq)).gt.tresh)then
  476.  
  477.               h=d(iq)-d(ip)
  478.  
  479.               if(abs(h)+g.eq.abs(h))then
  480.  
  481.                 t=a(ip,iq)/h
  482.  
  483.               else
  484.  
  485.                 theta=0.5*h/a(ip,iq)
  486.  
  487.                 t=1./(abs(theta)+sqrt(1.+theta**2))
  488.  
  489.                 if(theta.lt.0.)t=-t
  490.  
  491.               endif
  492.  
  493.               c=1./sqrt(1+t**2)
  494.  
  495.               s=t*c
  496.  
  497.               tau=s/(1.+c)
  498.  
  499.               h=t*a(ip,iq)
  500.  
  501.               z(ip)=z(ip)-h
  502.  
  503.               z(iq)=z(iq)+h
  504.  
  505.               d(ip)=d(ip)-h
  506.  
  507.               d(iq)=d(iq)+h
  508.  
  509.               a(ip,iq)=0.
  510.  
  511.               do 16 j=1,ip-1
  512.  
  513.                 g=a(j,ip)
  514.  
  515.                 h=a(j,iq)
  516.  
  517.                 a(j,ip)=g-s*(h+g*tau)
  518.  
  519.                 a(j,iq)=h+s*(g-h*tau)
  520.  
  521. 16            continue
  522.  
  523.               do 17 j=ip+1,iq-1
  524.  
  525.                 g=a(ip,j)
  526.  
  527.                 h=a(j,iq)
  528.  
  529.                 a(ip,j)=g-s*(h+g*tau)
  530.  
  531.                 a(j,iq)=h+s*(g-h*tau)
  532.  
  533. 17            continue
  534.  
  535.               do 18 j=iq+1,n
  536.  
  537.                 g=a(ip,j)
  538.  
  539.                 h=a(iq,j)
  540.  
  541.                 a(ip,j)=g-s*(h+g*tau)
  542.  
  543.                 a(iq,j)=h+s*(g-h*tau)
  544.  
  545. 18            continue
  546.  
  547.               do 19 j=1,n
  548.  
  549.                 g=v(j,ip)
  550.  
  551.                 h=v(j,iq)
  552.  
  553.                 v(j,ip)=g-s*(h+g*tau)
  554.  
  555.                 v(j,iq)=h+s*(g-h*tau)
  556.  
  557. 19            continue
  558.  
  559.               nrot=nrot+1
  560.  
  561.             endif
  562.  
  563. 21        continue
  564.  
  565. 22      continue
  566.  
  567.         do 23 ip=1,n
  568.  
  569.           b(ip)=b(ip)+z(ip)
  570.  
  571.           d(ip)=b(ip)
  572.  
  573.           z(ip)=0.
  574.  
  575. 23      continue
  576.  
  577. 24    continue
  578.  
  579.       pause 'too many iterations in jacobi'
  580.  
  581.       return
  582.  
  583.       END
  584.  
  585.     SUBROUTINE dcholdc(a,n,np,p)
  586.  
  587. C  Numerical Recipes Software 
  588.  
  589. c   Corrected by schuenemeyer 3/4/94.
  590.  
  591. c   Cholesky square root decomposition of n x n matrix.
  592.  
  593. c   Stored in lower triangular of a.
  594.  
  595. c   Double precision version
  596.  
  597.       INTEGER n,np
  598.  
  599.       REAL*8 a(np,np),p(n) ,sum
  600.  
  601.          
  602.  
  603.       do 13 i=1,n
  604.  
  605.         do 12 j=i,n
  606.  
  607.           sum=a(i,j)
  608.  
  609.         if(i.gt.1) then
  610.  
  611.           do k=i-1,1,-1
  612.  
  613.             sum=sum-a(i,k)*a(j,k)
  614.  
  615.           end do
  616.  
  617.         end if
  618.  
  619.           if(i.eq.j)then
  620.  
  621.             if(sum.le.0.) then
  622.  
  623.            write(*,*) 'choldc failed',i,sum
  624.  
  625.            pause 9
  626.  
  627.           end if
  628.  
  629.             p(i)=dsqrt(sum)
  630.  
  631.             a(i,i)=p(i)
  632.  
  633.           else
  634.  
  635.             a(j,i)=sum/p(i)
  636.  
  637.           endif
  638.  
  639. 12      continue
  640.  
  641. 13    continue
  642.  
  643.       return
  644.  
  645.       END
  646.  
  647.     SUBROUTINE BUBLE(X,ID,N)
  648.  
  649. c   Sorts X in ascending order, carries along ID
  650.  
  651. c   X is real, N is number of observations in X.
  652.  
  653.       DIMENSION X(1),ID(1)
  654.  
  655.       KS=N
  656.  
  657.    15 KW=0
  658.  
  659.       DO 30 I=2,KS
  660.  
  661.       IF(X(I).GE.X(I-1)) GOTO 30
  662.  
  663.       TMP=X(I)
  664.  
  665.       X(I)=X(I-1)
  666.  
  667.       X(I-1)=TMP
  668.  
  669.       NTI=ID(I)
  670.  
  671.       ID(I)=ID(I-1)
  672.  
  673.       ID(I-1)=NTI
  674.  
  675.       KW=1
  676.  
  677.    30 CONTINUE
  678.  
  679.       IF(KW.EQ.0) RETURN
  680.  
  681.       KS=KS-1
  682.  
  683.       IF(KS.EQ.1) RETURN
  684.  
  685.       GOTO 15
  686.  
  687.       END
  688.  
  689.     SUBROUTINE BUBLEI(X,ID,N)
  690.  
  691. c   Sorts X in ascending order, carries along ID
  692.  
  693. c   X is integer, N is number of observations in X.
  694.  
  695.       integer x,tmp
  696.  
  697.       DIMENSION X(1),ID(1)
  698.  
  699.       KS=N
  700.  
  701.    15 KW=0
  702.  
  703.       DO 30 I=2,KS
  704.  
  705.       IF(X(I).GE.X(I-1)) GOTO 30
  706.  
  707.       TMP=X(I)
  708.  
  709.       X(I)=X(I-1)
  710.  
  711.       X(I-1)=TMP
  712.  
  713.       NTI=ID(I)
  714.  
  715.       ID(I)=ID(I-1)
  716.  
  717.       ID(I-1)=NTI
  718.  
  719.       KW=1
  720.  
  721.    30 CONTINUE
  722.  
  723.       IF(KW.EQ.0) RETURN
  724.  
  725.       KS=KS-1
  726.  
  727.       IF(KS.EQ.1) RETURN
  728.  
  729.       GOTO 15
  730.  
  731.       END
  732.  
  733.